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Resume 

Pour un probleme elliptique bidimensionnel, nous proposons de formuler la 
methode des volumes finis avec des elements finis mixtes de Petrov-Galerkin qui 
s'appuient sur la construction d'une base duale de Raviart-Thomas. 

Abstract 

For an elliptic problem with two space dimensions, we propose to formu- 
late the finite volume method with the help of Petrov-Galerkin mixed finite ele- 
mentsthat are based on the building of a dual Raviart-Thomas basis. 
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1) Introduction. 

• Let fl be a bidimensional bounded convex domain in H 2 with a polygonal 
boundary dO,. We consider the homogeneous Dirichlet problem for the Laplace 
operator in the domain O, : 

(1.1) -Aw = / in Q 

(1.2) u = on the boundary dfl of Q. 

We suppose that the datum / belongs to the space L 2 (f2), 

(1.3) /el 2 (0), 

where this Hilbert space is classically defined according to 
L 2 (Q) = | v : Q — y H , / \v\ 2 dx < oo } 

(it, v) = / uv dx , Vu, v G L 2 (ri) . 

II w llj = (w, u) , Vm G L 2 (^). 

We introduce the momentum p defined by 
(1.5) p = Vu. 

Taking the divergence of both terms arising in equation (1.5), taking into account 
the relation (1.1) and the hypothesis (1.3), we observe that the divergence of 
momentum p belongs to the space L 2 (Q). For this reason, we introduce the 
vectorial Sobolev space H(div, Q) : 

#(div, O) = {qeL 2 (Q) xL 2 (Q), dwqeL 2 (Q)} 

H< (div>n) = /jkl 2 + |div,| 2 ]d,, V, eH(dW,n) 

and we suppose in the following that the momentum p satisfies the condition 
(1.7) p G i7(div, Q). 

• The variational formulation of the problem (1.1) (1-2) with the help of the 
pair £ = (u, p) is obtained by testing the definition (1.5) against a vector valued 



(1.6) 
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function q and integrating by parts. With the help of the boundary condition, it 
comes : 

(1.8) (p, q) + (u, dwq) = 0, G H(dw, Q) . 

Independently, the relations (1.1) and (1.5) are integrated on the domain Q after 
multiplying by a scalar valued function v G L 2 {Q) . We obtain : 

(1.9) (divp, v) + (/, v) = 0, Vv£L 2 (Q). 

The "mixed" variational formulation is obtained by introducing the product space 
V defined as 

f V = L 2 (Q) x #(div, ft), 
\ II 0, p) \\ v = || u \\ 2 + || p \\ 2 + || divp \\ 2 Q , 



(1.10) 



the following bilinear form 7(0, •) defined on V x V : 

(1.11) 7((w, p), (v, g)) = (p, g) + (w, divg) + (divp, v) 
and the linear form cr(«) defined on V according to : 

(1.12) <<r, C>= -(/,"), C = 0>,«) ev. 

Then the Dirichlet problem (1.1) (1.2) takes the form : 

(L13) 17(1,0 =<^,C>, v(ev. 

Due to classical inf-sup conditions introduced by Babuska [Ba71], the problem 

(1.13) admits a unique solution £ G V . 

• We introduce a mesh T that is a bidimensional cellular complex (see e.g. 
Godbillon [Go71]) composed in our case by triangular elements K (K G £ T ), 
straight edges a (a G A^) and ponctual nodes S (S G 5^). We conside also 
classical finite dimensional spaces L^ft) and ii^div, ft) that approximate the 

spaces L 2 (ft) and H(div, ft) respectively. A scalar valued function v G L^-(O) 
is constant in each triangle K of the mesh : 

(1.14) L 2 {SI) = {v.Q — > H, VK G £ r , 3v K G IR, Vir G if, u(a:) = v K ) ■ 

A vector valued function function q G ^^(div, ft) is a linear combination of 
Raviart- Thomas [RT77] basis functions ip a of lower degree, defined in the forth- 
coming section. 

• Let a G Aj- be an internal edge of the mesh, we denote by S and N the two 
vertices that compose its boundary da (see Figure 1) : 

(1.15) da = {S, N} 

and by K and L the two elements that compose its co-boundary d c a 

(1.16) d c a = {K, L] 
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in such a way that the normal direction n a is oriented from K towards L and 
that the pair of vectors (n a , is direct, as shown on Figure 1. We denote by 
W (respectively by E) the third vertex of the triangle K (respectively of the 
triangle L) : 

(1.17) K = (S, N, W) , L = (N, S, E) . 




Figure 1. Co-boundary (K, L) of the edge a = (S, N). 

The vector valued Raviart-Thomas [RT77] basis function ip a is defined by the 
relations 

1 [a: - W) , x eK 



;i-18) ipa(x) 



2 \K 
1 



2 \L\ 



(ar-E) 



I 



x G L 
elsewhere. 



When the edge a is on the boundary dfl , we suppose that the normal n points 
towards the exterior of the domain, so the element L is absent. We have in all 
cases the if (div, Q) conformity : 

(1.19) (p a G H(div, Q) 

and the degrees of freedom are the fluxes of vector field ip a for all the edges of the 
mesh (see [RT77]) : 

(1.20) / ip a »n a d7 = S aib , Va,beA T . 

J b 

A vector valued function q G i? T (div, Q) is a linear combination of the basis 
functions ip a : 

(1.21) q= Q*¥a e H T (div,n) =<ip b ,beA T >. 



• The mixed finite element method consists in choosing as discrete linear space 
the following product : 
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(1.22) V T = L 2 r {Q) x # r (div, Q) 

and proposes to replace the letter V by V T inside the variational formulation 



(1.13) : 

(L23) {tTc^cT = <<t,C>, V(gK 



T 

or in other terms 

( u T eL 2 r (Q) , p T e# r (div, Q) 
(1.24) i (p r , q) + (w r , divg) = 0, £ # r (div, ft) 

( (divp r ,^) + =0, V*; GL^fft). 

The numerical analysis of the relations between the continuous problem (1.13) and 
the discrete problem (1.23) as the mesh T is more and more refined is classical 
[RT77]. The above method is popular in the context of petroleum and nuclear 
industries but suffers from the fact that the associated linear system is quite dif- 
ficult to solve from a practical point of view. The introduction of supplementary 
Lagrange multipliers by Brezzi, Douglas and Marini [BDM85] allows a simplifica- 
tion of these algebraic aspects, their interpretation by Croisille in the context of 
box schemes [Cr2k] gives a good mathematical foundation of a popular numerical 
method and the possibility to reduce the size of the linear system has also been 
explored by Younes, Mose, Ackerer and Chavent [YMAC97]. 

• From a theoretical and practical point of view, the resolution of the system 

(1.24) can be conducted as follows. We introduce the mass-matrix M a ^ associ- 
ated with the Raviart-Thomas vector valued functions : 

(1.25) M a>b = (<p a , cp b ) , a,beA r . 

Then the first equation of (1.24) determines the momentum 

(1.26) p T = J2 PxaV* 

as a function of the mean values u TR for K £ Sj. : 



b eA T k e£ T 



K 



The representation (1.27) suffers at our opinion form a major defect : due to the 
fact that the matrix M _1 is full, the discrete gradient p T is a global function of 
the mean values u TR and this property contradicts the mathematical foundations 
of the derivation operator to be linear and local. An a posteriori correction of 
this defect has been proposed by Baranger, Maitre and Oudin [BM096] and with 
an appropriate numerical integration of the mass matrix M, it is possible to lump 
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it and the discrete gradient in the direction n a of the edge a is represented by a 
formula of the type : 

(1.28) p T = U%L - U%K 

V 1 ^ h a 

with the notations of Figure 1. The substitution of the relation (1.28) inside 
the second equation of the formulation (1.24) conducts to a variant of the so- 
called finite volume method. In an analogous manner, the family of finite volume 
schemes proposed by Herbin [He95] supposes a priori that the discrete gradient in 
the normal direction admits a representation of the form (1.28). Nevertheless, the 
engineer intuition is not correctly satisfied by a scheme such that (1.28). The finite 
difference UT > L ~ UT > K must be a priori to be a good (strong ?) approximation of 
the gradient p T = Vu T in the direction KL whereas the coefficient p Ta is an 

approximation of J a Vu T »n dr in the normal direction (see again the Figure 
1). When the mesh T is composed by general triangles, this approximation is 
not completely satisfactory at our opinion and contains a real limitation of these 
variants of the finite volume method. 

• We recall here that the finite volume method for the approximation of the dif- 
fusion operators has been first proposed from empirical considerations. Following 
e.g. Noh [No64] and Patankar [Pa80], the idea is to represent the normal interface 
gradient J a Vu T »n dr as a function of neighbouring values. Given an edge 
a, a vicinity V(a) is first determined in order to represent the normal gradient 
P Ta = J a Vu T »n dr with a "derivation formula" of the type 

(1.29) / Vu T .n dr = ^ g r r , u 



a K£V(a) 



a,K T,K 



Then the conservation equation 
(1.30) divp + / = 

is integrated inside each cell K G £ T is order to determine an equation for the 
mean values u TR for all K £ 8^-. The difficulties of such approches have been 

presented by Kershaw [Ke81] and a variant of such scheme has been first analysed 
by Coudiere, Vila and Villedieu [CVV99]. The key remark that we have done with 
F. Arnoux (see [Du89]), also observed by Faille, Gallouet and Herbin [FGH91] is 
that the representation (1.29) must be exact for linear functions u T . We took this 
remark as a starting point for our tridimensional finite volume scheme proposed 
in [Du92]. It is also an essential hypothesis for the result proposed by Coudiere, 
Vila and Villedieu. 



• In this contribution, we propose to discretize the variational problem (1.13) 
with the Petrov-Galerkin mixed finite element method, first introduced by Thomas 
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and Trujillo [TT99] . In the way we have proposed in [Du2k] , the idea is to construct 
a discrete functional space H T (div, 0) generated by vectorial functions </?*, a G 
A T , that are conforming in the space H(div, Q) 

(1.31) & e H(dw, n) 

and to represent exactly the dual basis of the family b G A T } with the L 2 
scalar product : 

(1.32) {ipa^t) = S a , b , Va, b £A T . 

(1.33) #*(div, fi)=<^,6GA r >C #(div, ft) . 

Then the mixed Petrov-Galerkin mixed finite element method consists just in re- 
placing the space i?,^ (div, ft) by the dual space H T (div, ft) for test functions 
in the first equation of discrete formulation (1.24). We obtain by doing this the 
so-called Petrov-Galerkin finite volume scheme : 

f u T GL^(ft) , p T G# r (div, ft) ^ 

(1.34) <^ (p T ,q) + {u T1 divq) = 0, V q G i?^- (div, O) 
[ (divp r ,^) + (/, V ) =0, Vt; GL^(fi). 

We introduce a compact form of the previous mixed Petrov-Galerkin formulation 
with the help of the product space V* defined by 

(1.35) V* = L 2 T (Q) x #*(div, ft). 
Then the formulation (1.34) admits the form : 

(L36) \7« T> =<^,C>, VCG7 r *. 

By doing this choice, it is easy to check that the scheme (1.34) is in fact a finite 
volume scheme for the Laplace operator. The key point is to construct the so-called 
dual Raviart-Thomas basis functions <p* in order to guaranty Babuska's 
[Ba71] inf-sup stability property. 

• The plan of the article is the following : we derive in the second part sufficient 
conditions in order to guaranty the final stability of the finite element scheme. 
Then we propose a particular family of dual Raviart-Thomas functions and propose 
by doing this a two-parameter family of finite volumes schemes. 



2) Stability analysis 

• We suppose in the following that the mesh T is a bidimensional cellular 
complex composed by triangles as proposed in the first section. Following the 
work of Ciarlet and Raviart [CR72], for any element K G Sj. we denote by h R 
the diameter of the triangle K and by p K the diameter of the inscripted ball 
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inside K. We suppose that the mesh T belongs to a family U Q of meshes that 
satisfies the following definition. 

Definition 1. Family of regular meshes 

Let 9 be a strictly positive parameter. The family U g of meshes is defined by the 
condition 

(2.1) r eu 6 yK e s T , < e. 

Pk ★ 

We suppose also that the dual space H T (div, Q) constructed by the conditions 
(1.31), (1.32), (1.33) satisfies the following hypothesis. 

Hypothesis 1. Interpolation operator H T (div, fi) — > (div, Q) . 
We suppose that the mesh T belongs to the family U Q of Definition 1 and that 
the dual basis ip* is constructed in such a way that there exists a linear map- 
ping if ^ (div, Q) 3 q i — > U q G if,^ (div, fi) and strictly positive constants 
A, B, D, E that only depends on the parameter 9 such that we have the following 
estimations : 

(2.2) A\\q\\l < (q,Uq) , V<? G # r (div, fi) 

(2.3) || II g || < 5 || q || , Vge # r (div, Q) 

(2.4) || div II g || Q < D || divg || , Vg G # r (div, Q) 

(2.5) (divg, divllg) > £ || divg \\ 2 Q , Vg G if r (div, Q) . 
Proposition 1. Divergence lifting of scalar fields 

Let 9 be a strictly positive parameter. We suppose that the dual Raviart-Thomas 
basis satisfies the Hypothesis 1. Then there exists some strictly positive constant 
F that only depends on the parameter 9 such that for any mesh T that belongs 
to the family U Q , and for any scalar field u constant in each element K G £j. 
(u G L^_(Q)), there exists some vector field q G if ^ (div, Q) such that 

( 2 - 6 ) Il9ll*(div,n) ^ F H^llo 

(2.7) («,divg) > \\u\\l . 
Proof of Proposition 1. 

• Let u G L^_(fi) be a discrete scalar function supposed to be constant in each 

triangle K of the mesh T. Let i/j G Hq(Q) be the variational solution of the 
Poisson problem 

(2.8) Ai^ = n infi, ^ = on^. 

Since Q is convex, the solution i/j of the problem (2.8) belongs to the space if 2 (f2) 
and there exists some constant G > that only depends on Q such that 
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(2.9) || if) || 2 < G || u || . 

• Then the field belongs to the space H X (Q) x H 1 (Q). It is in consequence 
possible to interpolate this field in a continuous way (see e.g. Roberts and Thomas 
[RT91]) in the space H(div, Q) with the help of the fluxes on the edges : 

(2.10) Pa = I ^d 7 , P = P a V a ^ H T (div,Q) 

and there exists a constant L > that only depends on the parameter 9 such 
that 

( 2 - n ) IMItf(div,o) < L \\u\\ Q . 

• We observe that we have exactly 

(2.12) divp = u in O. 

On one hand, the two fields divp and u are constant in each element K of the 
mesh T. On the other hand, we have : 

[ divp dx = f p»n d 7 = /" dj = f Atp dx = / u da: 

J K JdK JdK V n J K J K 

and the relation (2.12) is a consequence of the above property for the mean values. 

• Let lip be defined according to the Hypothesis 1, and 

(2.13) q = in P . 

We have as a consequence of (2.5) and (2.12) : 

(u,divq) = — (divp, divLTp) > || divp |P = || u \\ 2 
E oo 

that establishes (2.7). Moreover, we have due to (2.3), (2.4) and (2.11) : 

1 II XT M B M M BL 

M y no H "o E ~~ 



5 




E 


IIP Mo ^ 




D 







~ E ll 



i , d . , d „ 

||divg|| = - ||divnp|| < — ||divp|| = — \\u\\ . 

Then due to the definition (1.6), the two above inequalities establish the estimate 
(2.6) with F = -gr y/B 2 L 2 + D 2 and the Proposition is proven. □ 

Proposition 2. Discrete stability 

Let # be a strictly positive parameter. We suppose that the dual Raviart-Thomas 
basis satisfies the Hypothesis 1. Then we have the following discrete stability for 
the Petrov-Galerkin mixed formulation (1.36) : 

f 3/3> 0, VTe U e , V£ eV T such that H\\ v =l, 

[ } I 3r/ € V£, || || v < 1 and 7 (£, »/) > • 



FRANgois Dubois 



with 7(», •) defined at the relation (1.11) and (3 chosen such that 

(-5) f^W 7 ^ >- ( i+F ( i+ ^¥))^- 

Proof of Proposition 2. 

• We set £ = (w, p) satisfying the hypothesis (2.14) : 

(2.16) n\\ 2 v E= huh; + hpii; + iidivpyj = 1. 

Then at last one of these terms is not too small and due to the three terms that 
arise in relation (1.11), the proof is divided into three parts. 

• If the following condition 

(2.17) ||divp|| > /3, 
is satisfied, we set 

(2 - 18) " = tSt ' q = ' i = ^^ e K- 

Then || div-y || Q = 1 and || £ || < 1. Moreover 7(£, C) — (divp,v) = 
= || divp || > (3 by hypothesis (2.17) and the relation (2.14) is satisfied in this 
particular case. 

• Under the following conditions 

(2.19) ||divp|| < f3 and || p \\ Q > \/ - y/p , 



we set 



(2.20) v = , q = -^-^ Up, ( = («, q) € V* . 

The following inequalities are a direct consequence of (2.3) and (2.4) : 

b .. „ .. d 

Q L < ^ P L an d divo L < — — divp L 

ii y " " o - b + D 11110 " "° _ B + D " "° 

so we deduce : 

ll2 B 2 + D 2 „ Il2 

because B > 0, .D > 0. Then || £ || < 1 and we have also 
7(£, C) = (P, 9) + (w, divg) + (divp, v) 

^ 1 1 „2 ^ n „ 

> ^ D ((B + 2D)P - D(3) because || u \\ Q < 1 

> (3 and the property is established in this case. 
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• If the two previous conditions (2.17) and (2.19) are in defect, i.e. if we have 

(2.21) || divp || < and || p \\ Q < ]/^^^ V? , 

B -\- ID 

then \\u \\l = 1 - II p \\l - || divp |P > 1 /3 - f3 2 > f3 > 

ii iiq ii x- iiq ii 11 A i — i 

due to £he hypothesis (2.15). From the Proposition 1, there exists some vector field 
q E H T (div, fi) satisfying || q || H(diVj Q) < || w || and (w,divg) > 



u\\ < 1 

due to the hypothesis (2.14) relative to || £ || y . Moreover, we have 
7(f, C) = (P, 9) + divg) + (divp, v) 

= -p (P, 9) + -p (w, divg) 



ii - 112 


We set 




(2.22) 


v = , q = 


1 

F 


then 








1 




IICIIv 


— p II Q 1 1^ (div.n) 


< 



If B + 2D r- „ „ „ 




/ 1 „ B + 2D r- 

> IMIo H^llo ~V A ^ 

> /3 and the property is satisfied for this last case. 
The Proposition 2 is established. □ 

Theorem 1. Error estimate 

Let f2 be a two-dimensional open convex domain of IR 2 with a polygonal bound- 
ary, u E H 2 (Q) be the solution of the problem (1.1) (1.2) considered under varia- 
tional formulation and p = Vu be the associated momentum. Let 9 be a strictly 
positive parameter, U Q a family of meshes T and V* defined in (1.35) and asso- 
ciated with a choice of a dual Raviart-Thomas basis that satisfies the Hypothesis 
1. Let £ = {u T , p T ) E Vj. be the solution of the discrete problem (1.34). Then 
there exists some constant C > that only depends on the parameter 6 such 
that 
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(2.23) ||«-« r || + \\p-p T ll H(div>n) < Ch T ||/|| . 
Proof of Theorem 1. 

• On one hand, it is sufficient to apply the general approximation Theorem 
established by Babuska's for continuous (respectively discrete) variational mixed 
systems (1.13) (respectively (1.36)) i.e. to verify that the bilinear form 7 («, •) 
defined in (1.11) is continuous on the Hilbert space V = L 2 (Q) x H(div, fi) , 
which is clear. It is also necessary to verify the so-called discrete inf-sup condition 
(2.14), that has been established at the Proposition 2. Last but not least, it is 
necessary to satisfy the following infinity condition : 

(2.24) Vr/EV* sup 7 (£, 77) = +00 . 

1 dev T 

• The infinity condition (2.24) is established as follows. Let £ = (v, q) E be 
a "test vector" different from zero. If there exists some mesh element K E ^ such 
that J^divgda; 7^ 0, then we consider £ = (it, p) chosen according to u = Xu 
and p = 0. We suppose that the field it E is null for all the elements of the 

mesh T except for the particular element K where we suppose u R = J K divqdx. 
Then we have 7 (£, £) = (p, + divg) + (divp, v) = A (J^divgdaf) , which 
tends to infinity as A tends to infinity. If J K div g da: = for all mesh elements 
K E Ej. and if the field q is not null, we can write it on the form q = Up with 
p E if^div, Q) because the mapping II is clearly bijective due to the property 
(2.2). We set p = Ap and it = 0. Then 7 (£,C) = (p, q) = A (p , q) > XA \\ 
q 11^ due to the hypothesis (2.2) ; the infinity property (2.24) is established in this 
second particular case because q ^ 0. If q = 0, then v is not null due to the left 
hand side of (2.24). Following the proof of Proposition 1, we introduce the vector 
field p E .f/^div, fi) satisfying the relations (2.11) and (2.12) : divp = v and 
|| p ^ < L \\ v || . We set p = Ap, it = and £ = (it, p). Then 

7(£>C) = A (divp, i>) = A || 1; 11^ tends to infinity when A tends to infinity, and 
the infinity condition (2.24) is established. 

• The conclusion of the Babuska's Theorem [Ba71] assures the existence of some 
constant C > that only depends on 6 such that the error between the solution 
of the continuous problem (1.13) and the discrete problem (1.23) is maj orated by 
the interpolation error : 

r \\u-u T n + 11 p-p t n H(div>n) < 

(2.25) < < C ( inf || it - v L + inf \\p-q\\rr (A - J. 
( ° <?Etf r (div, fi)" "^'^ 

Then following classical interpolation results for scalar [CR72] and vectorial [RT77] 
fields, we deduce from (2.25) : 
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II u-u T || + \\p-p r || H(diV)0) < C (h T \\u\\ x +h T HpIIJ 

< C /i r || w || 2 < C h T || / || 
and the Theorem 1 is proven. □ 



3) Towards a first Petrov-Galerkin finite volume scheme 

• We propose in this section to formulate some ideas in order to construct a dual 
Raviart- Thomas basis tp* where a is an internal edge of the mesh T {a £ A T ). 
Following (1.15) and (1.16), we denote by a = (S, N) this edge, by O the middle 
of SN and by K, L the two triangles that compose the co-boundary. The normal 
n a is supposed to be oriented from the element K towards the element L and 
there exists two vertices W and E such that the relation (1.17) holds (see the 
Figure 1). We consider the four edges (N, W), (W, S), (S, E) and (E, N) that 
compose the boundary of the union KUL. We define four new triangles M, P, Q 
and R and four new vertices A, B, C and D in the mesh T by the relations 





f 6> C (E, N) = (L, M) , 


M = 


(N, E, A) 


(3.1) < 


d c (N, W) = (K, P) , 


P = 


(W, N, B) 


d c (W, S) = (K, Q) , 


Q = 


(S, W, C) 




{ d c (S, E) = (L, R) , 


R = 


(E, S, D) 


as illustrated on the Figure 2. 








D 



Figure 2 : support V(S, N) of the dual Raviart-Thomas basis function ip* 
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Hypothesis 2. A simple choice for dual Raviart Thomas basis func- 
tions 

We suppose first that the Raviart Thomas dual basis cp* (b G A T ) satisfies the 

H(d\v) conformity property (1.31) and the orthogonality (1.32). Moreover, we 
suppose that for each internal edge a = (S, N), the support of the dual Raviart- 
Thomas basis function <p* is included in a vicinity V(a) = V(S, N) composed 
by the six triangles K, L, M, P, Q and R introduced previously (see the Figure 
2) : 

J V(S, N) c KULUMUPUQUR, 
^• 2 ) { supp(<4J c V(S,N). 

We suppose also that the divergence field div ip* is constant in each triangle 
of the mesh : 

(3.3) divcp* G L; 2 (Q) , Va G A.j- . 




Figure 3 : Notations and orientations. 



Theorem 2. Necessary condition for a dual Raviart- Thomas basis. 

Let </?* N be a dual Raviart Thomas basis function satisfying the Hypothesis 2. 
We introduce the following fluxes accross the internal edges SN, EN, NW, WS 
and SE respectively : 
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(3.4) 



7] 



a 



7 



'SN 



^SN #n SN 



'EN 



'WS 



^SN #n EN 



^SN #n WS 



d7, 
d7, 
d7, 



5 



'NW 



^SN # n NW 



d7, 



'SE 



^SN* n SE 



d7 



Then we have the necessary conditions : 

(3.5) rjKt + aufr + + 7^ + 5Ul = |Sft| 



n 



SN 



(3.6) 



a 



UJl.WA + /3KP\EF3 + ~/K$.E& + <5LFi.Wr3 = 

= -3 |SN| n SN .(ot + Ojt) 



• The finite volume approach is then obained with a six point scheme for the 
mean gradient in the normal direction in the manner of (1.29) thanks to the first 
equation of the mixed variational formulation (1.24) : 



(3.7) ( / SN V ^ #n d7 = ^l-^kH 
[ +a(u M -u h ) + f3(u p -u K ) - 



+ 7(w Q - u K ) + 5(u R - u L ). 

We remark that the constraints (3.5) express that the relation (3.7) is exact if the 
field Uj. is an affine function. 




Figure 4 : Notations for an arbitrary edge BE. 

• We precise some notations that we will use in the next pages. Let (B, E) 
be an edge of the mesh (see e.g. the Figure 4), and (L, R) its co-boundary. If 
the edge is directed from B towards E, the axis s has its origin at vertex B 
and the normal n BE is oriented from L to R in such a way that the pair of 
vectors (n BE , BE) is direct. If £ = J BE tp* »^ BE ds is the flux^of the function 
cp* accross the edge (B, E), we will denote by £, , £ 1 , £ 2 , and £ 2 the following 
momenta : 
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(3.8) 



fi = / l P** n BE sds ' ^ = ^** n BE (BE-s)ds, 
Jbe Jbe 



'BE 



(p*.n BE s 2 ds, f. 



'BE 



<p*»n BE (BE-s) 2 ds 



The proof of Theorem 2 needs a certain number of technical lemmae and prelimi- 
nary propositions. 

Lemma 1. Radius of giration 

Let M = (N, E, A) be a triangle of the mesh T and M its associated center of 
gravity (see the Figure 4). We will denote by p M the radius of giration : 

(3.9) p M = y i (AN 2 + NE 2 + EA 2 ) 
Then we have 

(3-10) f M k-N| 2 dx = p 2 M +NM 2 . 




x 



Figure 5 : About the radius of giration. 



Proof of Lemma 1. 

We have on one hand 

\M\= ^h{p + q), 
/ \x - N| 2 dx = 

J M 



■h 



'0 



dy 



y 

[ h y dx (x 2 + y 2 ) 
J—v — 



■ h 



(p + q) {P 2 - PQ + Q 2 ) h + | {p + q) h 3 
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1 1 

= \M\ [-{p 2 -pq + q 2 ) + -h 2 ) 
and on the other hand : 
i (AN 2 + NE 2 + EA 2 ) + NM 2 = 

= ^[p 2 + h 2 +q 2 + h 2 + {p + q) 2 } + \[{p-q) 2 + ±h 2 } 

= h ( 6 p 2 + 6 ^ 2 - 6pq + 18ft 2 ) 

36 

= \ {p 2 + q 2 - pq + 3/i 2 ). 

So the relation (3.10) is established. □ 



Proposition 3. First relations between momenta 

The Hypothesis 2 implies the following relations inside the triangle M = (N, E, A) 



^ri Eft Eft 

ft =EM. — a, a, = -NM. — a, 

1 EN ' 1 EN ' 

<*2 = (Pi + EM2 )« > «2 = (PL + NM2 )« 



(3.11) 

and the analogous ones obtained from the Figure 3 in the triangles P = (W, N, B), 
Q = (S, W, C) and R = (E, S, D) : 




— * N# ~ — i NV^ 

= NP . B , & = -WP . /3 , 

NW ^ ' ^ NW M ' 

= (Pp +NP 2 )/3, /3 2 = ( P 2 p + WP 2 )/3, 



WS -4 ws 

71 = ^'ws 7 ' 7l = _sq, ws 7 ' 

7 2 = (P 2 Q + WQ 2 ) 7 , 7 2 = (P 2 +SQ 2 ) 7 , 



(3-12) 



(3.13) 



\ 6, = SR. — 5 , £ = -ER.— (5, 

(3.14) <^ 1 SE J SE 

1 S 2 = (p 2 R + SR 2 )<5, <5 2 = (p 2 + ER 2 )<5. 

Proof of Proposition 3. 

• We write the orthogonality (1.32) between the two edges a = (S, N) and the 
edge b = (A, N) (see the Figure 3). Inside the triangle M = (N, E, A), we have 

^an = - T^ V(|a; - E|2) 

then 
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= L * ^AN dx = / ^SN ' ^AN d * 

JQ, J M 



^^IW\ lx - E{2dx + L { ^' n) TW\ |a; " E|2d7 



JM t|JKI| JQ M 



1 




4 |M 



K N #n N E ) l^- E | 2 d 7 



NE 

£-E| 2 d:r- / (<p* N .n NE ) \x-E\ 2 ds 



M\ J M 1 1 Jen ^ 



a (p 2 M + EM 2 ) - a 2 ) 



4 |M 

and the third relation of (3.11) is proven. 

• We write now the orthogonality (1.32) between the two edges a = (S, N) and 
b = (E, A) inside the triangle M = (N, E, A). When we exchange the roles of the 
two vertices N and E in the previous relations, we obtain the same result, excepts 
that a 2 has to be replaced by a 2 . So the fourth relation of (3.11) is established. 

• We have from the relation (3.8) : a 2 = NE 2 a — 2 NE a 1 + a 2 . Then 

( EM 2 - NM 2 + NE 2 ) Eft . ( EM" + NM" + Eft ) 

ai = ZZ _v 2NE a = 2NE a 

EN. EM 

= a 

NE 

and the first relation of (3.11) is established. As previously, the exchange of the 
vertices N and E induces the change of a 1 into a 1 that establishes the second 
relation of (3.11). 

• The relations (3.12), (3.13) and (3.14) are obtained by circular permutation, 
following the rules that are natural when viewing the Figure 3 : 

a — ► 13 — ^7 — >5, E^N^W — >S, N — ^W — ^S^E 

and M — ^ P — ^ Q — > R. □ 

Lemma 2. A mean value of the dual Raviart- Thomas basis function 

Let M = (N, E, A) be a triangle of the mesh T associated to the edge a = (S, N) 
as in Figure 3 and </?* N a dual Raviart- Thomas basis function satisfying the 
Hypothesis 2. Then for each constant vector £, we have : 

(3.15) £. f <p* dx = (^.n EN )(M^.n EN )a. 

JM 
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Proof of Lemma 2. 

We have : 



J M J M 

= - [ ^.(x-M)div^* 6x + f (p* .n)£.(s-M)d7 

JM JdM 

= -(div^)/ £.0r-M)dr + / (<p* .(-n EN ))^.(^-M)d 7 



+ 



/ K N -(-^E N ))^b- E + E - M ] d 7 

JEN 



/ (v 

JEN 



SN #n EN' C * 



Ert 
EN 



+ m4 



d7 



V EN 1 / s [ V EN/ EN 

= (M^.n EN ) n EN a 

and the relation (3.15) is established. 



□ 



Lemma 3. A simple relation between two triangles 

Let L = (S, E, N) and M = (N, E, A) be the two triangles of the mesh T 
associated to the edge a = (S, N) as in Figure 3. Let ip* be the dual Raviart- 
Thomas basis function satisfying the Hypothesis 2. Then we have : 



( 3 - 16 ) 2\M\ L {A ~ X)9<P * n6X = 2|T| ( ^ #n EN)^- n E N) 



a 



Proof of Lemma 3. 

We denote by h the height of the triangle NEA and by h* the height of ENS 
chosen such that | M | = | h NE and \L\ = \h* NE. We have from the Lemma 
2 : 



2| 



1 -m\L {a - x) -* 



* dr 
SN aX 



— r- — - • / 

2 M 7m 



SN 



(f c „ dx 



i-/.(Elt.n EN ) 



/iNE 



a 



/i* NE 



h* (Elt.n EN ) a 



and the relation (3.16) is proven 



a 



□ 
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Proposition 4. Second relation between momenta 

The Hypothesis 2 inside the triangles L = (S, E, N) and M = (N, E, A) implies 
the following relation : 

r o| + SL 2 )r/ + isl.(s^ + si + srt) a + 
(3 - 17) = [ + i^.(srUsE>s^. 

We have also : 1 

[ (p| + NL 2 )r/ + -s1.(n3 + NE" + Nl) a + 

(3 - 18) ^ = 6 + 1^.(5*3 + ^ + ^)* , 

1 6 

f (pl: + NK 2 )r/ - -SB.(NB + N# + N^)/3 

(3 - 19) ^ = ° Ii^.(N^ + Nd + N3) 7 , 

1 

[ (p^ + SK 2 )?? - -SB.(Srt + SB + SV^)/3 

(3 - 20) = { - + + 

Proof of Proposition 4. 

• We write the orthogonality (1.32) between the two edges a = (S, N) and the 
edge b = (E, N) (see the Figure 3). We have 

^-L-Or-S) = ^- V(|^-S| 2 ) inside L = (E, N, S) 

^EN = \ 1 

(A - a;) inside M = (N, E, A) . 



2 |M 



Then 

= / ^M-^N d ^ = / ^M^EN^ + / l P^* ( P^ dX 

J M 



. SN ^EN / ~SN ^ EN 1 / ^SN ^EN 



4 



J^(di^* N )(i)./ L l-^-SI 2 dx + ^L-J J^.n) \x-Sf d 7 



1 

2 \M\ J M 



+ ^FT / ¥>* N • (A - a) dr 



1 { 5 -\- a — rj 



4 L V LI 



a; - S| 2 da; + 
E^ 



+ ^TZT (g ^ -n EN)(^* n EN) « 
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and we have, thanks to the relations (3.10), (3.11) and (3.14) 



-{5 + a - 77) (p 2 L + SI/) + (p 2 R + SR 2 ) 5 + SE z a + 

(3.21) V 2 =<j +2 ( S ^.g)(ESi.g)a + (^ + EM> + 

+ 2(Srl.n EN )(E5t.» EN ) a. 

• The a coefficient in the right hand side of the relation (3.21) is equal to : 
-i(SE 2 + EN 2 + NS 2 ) - SL 2 + SE 2 + 2 (s4.Efti-(s4 .n EN ) 0t .n EN )) + 

+ i(EA 2 +AN 2 + NE 2 ) + EM 2 + 2 (SI&.n EN ) 0L .n EN ) 

= i ((El + sS).(El-Si) + (Aft + N^).(A^-N^)) - SL 2 + (S^ + Elt) 2 
36 

^s1.(e! + eS + n^ + n^) + (S^t + st) • (S^t - st) 

= — sl.(E^ + E^ + Nl + N^) + (2 St + Lit). Lit 
36 

= -Ls1.(e^ + e3 + n1 + n3) + (|(si + sify + isi).(isi) 

= ^sl.(El + E^ + Ni + N§ + 8(s4 + S^)+4Sl) 

= ^s1«(eS + s1 + eS + nS + s1 + nS + 8s4 + 8Srt + 4Sl) 

36 
6 

in coherence with the right hand side of the relation (3.17). 

• In a similar way, the coefficient of 5 in the right hand side of the relation 
(3.21) is equal to : 

-i(SE 2 + EN 2 + NS 2 ) + i(DE 2 + ES 2 + SD 2 ) + (sS + SL* ) . (Sit - st ) 
36 36 

= — ( (L^ + E^) . (L^ - E^) + (Sl5 + N^) • (Sl5 - N^) ) + 

ou 

+ i(SL5 + 2S^ + Srt).(iNL5) 

= — N6.(El5 + Elt + Sl1 + Srt + 4SL5 + 8Sfi + 4Srt) 
36 

= iN^.(E^ + ^ + E^ + ^ + Sr5 + Srt + 4Sl1 + 8S^ + 4S^) 
36 

= irrt.(sr5 + si + s^) 
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as proposed in the right hand side of the relation (3.17). Then the relation (3.17) 
is a direct consequence of (3.21). 

• The proof of the relation (3.18) is obtained from the previous relation (3.17) 
with the following changes : E i — ► E , rj i — > rj , A < — > D , N < — > S , M < — > R , 
a i — > 5 and r\ 2 < — > rj 2 . In a similar way, the relations (3.19) and (3.20) are 
a straightforward consequence of the relations (3.17) and (3.18) with a vision of 
the Figure 3 "from the top to the bottom" , id est with the following changes : 
E^ — ^W,N^ — ^S,D^ — ^B,A^C, L^K,M^Q,R^P, 
7] < — > —7] , a i — > 7 , 5 i — > (5 and r] 2 i — > —rj 2 . So the proposition is established. 
□ 

Proposition 5. Two expressions for the first order momentum 

Under the Hypothesis 2 and the notationa^Droposed at the Figure 3, we have : 

(3.22) j] x = (S£r/ + Llta + Llh).^- , 

(3.23) Vl = (S^rj-^P-K^^M . 



Proof of Proposition 5 . 

• We deduce from (3.8), (3.17) and (3.18) : 




and the relation (3.22) is established. 
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• The proof of the relation (3.23) is analogous. It is a consequence of the 
relations (3.8), (3.19) and (3.20) : 

1 

2NS 



v-i = ^Tatc ( % - % + sn2 "n ) 



( (SK 2 -NK 2 + SN 2 )r/ + \ SB . + + W§) f3 + 
z IN b V b 



+ ^NC5.(N^ + W^ + C^)7 + isF3.(NL^ + N# + N^)/3 + 
6 6 



+ ^NC5.(N# + NC5 + nS) 7) 

6 / 



2NS 



(Sit + Nit) . (Sit - Nit) + SN 2 ] rj + 

+ -^.(N^ + L^ + W^ + N^ + rTV^ + N^)/3 + 
6 

iN^.(N^ + W^ + C^ + N# + Nc5 + N^)7 



+ 



6 



(sN.(Spt + Nlt + SN) 77 + ?SB.N^/3 + ^Nc5.N^ 7 



2NS 



= ^(Sl4 77 - KP"/3 - K^ 7 )-SN. 

The relation (3.23) is established and the Proposition 5 is proven. □ 
Lemma 4. Two usefull integrals 

Let K = (S, N, W) and L = (N, S, E) be the two triangles of the mesh T that 
compose the co-boundary of the edge a = (S, N) as in Figure 3. Let </?* N be the 
a Raviart-Thomas basis function satisfying the Hypothesis 2. Then we have : 



(3.24) jL- I p* N . (E - z) dr = ^(s£ rj + uk a + LlU)... SN 



L 



(3 - 25) 2\K\j K ^ m{W ~ x)diX = ik^ 71 ~ ^ P ~ ^^)* n SN- 
Proof of Lemma 4. 

• We establish the relation (3.24) by integrating by parts and using the relations 
(3.10), (3.11) and (3.14) : 



2 



}l\ /^SN-( E -*) d * = -4TLI / L ^SN- V d*- E | 2 ) d * 

^(div^ SN ) ( L) J l |,-E| 2 ds - ^ Jj^n) |,-E| 2 d 7 
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+ 



+ 



i^.(SFi + A^ + sfi + srt + A^ + srt + 6(SFi + sl + srt)) 

36 

+ L]\t • (LFj + Mill + 2 Els) + 2 (SE* • n gN ) (L]\t • n gN ) a + 
iNl^.(N^ + L^ + N^ + ^ + D^ + SF$ + 6(Sr5 + SF$ + S^)) 

DO 



4 \L 



+ LFt.(L^ + + 2E^) + 2(S^.n SN )(Ltt.n SN ) 
1 /r 2(S^.n SN )(Slt.n SN )] 77 + 



S 



+ [i^.(2^+^+2^) + i^.(l^+M^) + 2(^.n SN )(LIvt.n SN )] a + 

1 1 

+ [-ND\(2S^ + ST^ + S^) + -ND\(l3 + r3) + 2(SE\n )(LFl.n SN )] 5 
9 3 " 



2 \L 
1 



(s4.n SN ) (s£.n gN )r/ + (ui*n SN ) a + (LFl.n SN ) 6 
(s£«n qN ) 77 + (Ot«n SN ) a + (LTl«n qN ) 5 



that establishes the relation (3.24). 

• The relation (3.25) is a consequence of the previous relation (3.24) with the 
following modifications : E i — > W , N < — > S , D < — > B , A < — > C , L < — > K , 
M i — > Q , R i — > P , 77 i — y -77 , a < — > 7 , 5 < — > (3 and n q < — > -n 



J SN 



□ 



Proof of Theorem 2. 

• We first eliminate the variable 7] 1 between the relations (3.22) and (3.23) ; 
we obtain 

(3.26) (Ktrj + Ilka + KP* /3 + KC^ 7 + Lit<5).Srt = 0. 

• We write secondly the ort honor mality relation (1.32) between the vector func- 
tion cp QN and its dual <£g N ? with the help of (3.24) and (3.25) : 

1 = (^SN'^SN) = / ^SN'^SN^ + / ^SN #( ^SN d * 



2 I K I J K ^sn v > 2 \L\ J L s N ' 

= -L ( (st + KlJ) rj + LSfca + KP^ + K^ + Llt 6). 

IN o 



n 



SN 
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J- ((st + Krt + NS) 7] + Lit a + KP*/3 + K^ 7 + Hl8).n 



NS 



SN 



= ^(K^ + + + K^ 7 + LlU).n SN . 

The relation (3.5) is a direct consequence of the above expression and of the 
previous relation (3.26). 

• Thirdly, we eliminate the variable r] 2 between the relations (3.17) and (3.20). 
The coefficient of the variable 77 is equal to 

(pi + SL 2 ) - (p 2 K + SK 2 ) = -L ((EN 2 + ES 2 ) - (NW 2 + SW 2 )) + SL 2 - SK 2 

36 

= i ((EfJ + N#).(EfJ + WfJ) + (E^ + S#).(E^ + W^)) + (SL + Slt).(SL + K^) 
= ^W^.(N^ + N# + S^ + S#) + Jw^.f J(S^ + S^ + S# + sS)) 

= -^-KL.(N^ + N# + SE" + S^ + 8Srt + 4SFJ + 4SM^) 

= -^-Kt.(N^ + sfi + N^ + S^ + SR* + S# + 8Srt + 4SE" + 4S^) 

J. ^ 

= ^KL.(S^ + Si + S#). 

We deduce from (3.17), (3.20) and the previous calculus : 

I^.(^ + ^ + ^)r/ + l^.(^ + ^ + ^)a + I^.(^ + SF3 + S#)^ + 

+ ^.(^ + S^ + S^) 7 + Llt.(Sl1 + si + S^) ( 5 = 

and taking into consideration the relation (3.26) : 

KL\ (SLJ + S\^) r/ + Lit . (sS + si) a + • (SF3 + S^) /3 + 

+ K^.(S# + S^) 7 + Lit. (sfe + si) £ = 0. 

We eliminate the variable r] be twee n the previous relation and the relation (3.5) 
after multiplying it by — (SE + sW). We get : 

Llvt.(si + Si-Si-S#)a + K^.(Sl1 + S#-Si-S^)/3 + 

+ K^.(S# + S^-s4-S^)7 + Lrt.(s6 + Si-Si-S^)5 = 

= -|NS| n gN .(si + S#) 

id est 

ui.WAa + KF" «EF3 ft + K^.E^ 7 + lS.WTj 6 = 

= -3 I NS I n gN . (SL + Sit) 

and the relation (3.6) is a direct consequence of the above relation. The Theorem 
2 is established. □ 
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4) Perspectives 

• We have proposed to formulate the finite volume method for the Poisson 
equation in two space dimensions with the help of Petrov-Galerkin mixed finite 
elements. The unknown is constant in each triangle and the momentum is dis- 
cretized with the Raviart- Thomas vectorial finite elements of lower degree. The 
conservation law is integrated in each triangle and our stencil for the discrete gra- 
dient operator is composed by six triangles in the vicinity of each edge of the mesh. 
The question of the determination of such a scheme conducts to a two-parameter 
family for a possible choice of a so-called "dual Raviart- Thomas basis function" 
for the finite volume scheme. We have also developed a sufficient hypothesis to 
prove the stability and the optimal convergence of the associated finite volume 
scheme. The next step of this research is to construct explicitly an interpolation 
vector valued function in the particular case where O, = JR 2 in order to determine 
free coefficients and to establish the stability property. 

• We thank Jean-Pierre Croisille for his kind invitation to first present [Du99] the 
results contained in this article and for regular helpfull discussions that convinced 
us of the complexity of mathematical links between the present formulation of the 
finite volume method and his analysis [Cr2k] of the box scheme. 
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